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Abstract 

In this paper we study the thermal effective behaviour for 3D multiphase composite material 
consisting of three isotropic phases which are the matrix, the inclusions and the coating media. 
For this purpose we use an accelerated FFT-based scheme initially proposed in Eyre and Milton 
(1999) to evaluate the thermal conductivity tensor. Matrix and spherical inclusions media are 
polymers with similar properties whereas the coating medium is metallic hence better conducting. 
Thus, the contrast between the coating and the others media is very large. For our study, we 
use RVEs (Representative volume elements) generated by RSA (Random Sequential Adsorption) 
method developed in our previous works, then, we compute effective thermal properties using 
an FFT-based homogenization technique validated by comparison with the direct finite elements 
method. We study the thermal behaviour of the 3D-multiphase composite material and we show 
what features should be taken into account to make the computational approach efficient. 
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I. INTRODUCTION AND MOTIVATIONS 


In this paper, we study the thermal conductivity tensor using the 3D homogenization 
method for materials constituted of 3 phases. Since composite materials are widely used 
in modern industry, there has been a number of works studying their effective properties: 
they concern descriptions of modelling approaches as well as results of experimental mea¬ 
surements. The purpose of this paper is to explore the range of application and features 
for existing modelling approaches, and at the same time the work is motivated by precise 
application within the framework of an industrial project. 

The importance of modelling for the analysis of composite materials results from various 
challenges in carrying out experimental works: such problems include purely technical dif- 
hculties as well as financial optimization of the applied activities. It is thus necessary to 
master modelling approaches that are reliable, comparable with experiment, and still af¬ 
fordable from the point of view of computational efficiency. Our strategy in this paper is 
related to the idea of stochastic homogenization. We consider samples of a composite ma¬ 
terial containing coated inclusions; and in order to take into account possible imperfections 
and random factors, we average the result for a series of tests representing the same macro 
characteristics. For this work we have chosen to study relatively simple geometry of inclu¬ 
sions: spherical shape with uniform coating. On the one hand this permits to concentrate 
on the details of application of the adopted modelling technique, on the other hand this 
can still provide some insight on the influence of morphology on the effective properties of 
composite materials. 

The paper is organized as follows. In the next section, we explain the RSA (Random Se¬ 
quential Adsorption) method chosen for the generation of RVEs (Representative volume 
elements) proposed in [T] and we give the context of the study. Next, we recall the model 
used with the computational method. In the third section we compare the FFT (Fast Fourier 
Transform)-based homogenization technique and the flnite elements method (FEM) for sim¬ 
ple geometries and we motivate our choice of the FFT-based ‘accelerated scheme’. The last 
section describes the numerical results. We conclude by recapitulating the main points of 
this work and describing the work in progress. 
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II. SAMPLE GENERATION AND SETTING OF THE STUDY 


In this section, we describe the method used to generate RVEs with three phases. As 
we have outlined in the introduction, we consider a material with the matrix, the spherical 
inclusions and the coating media. 

We proceed in two steps. First, we generate RVEs using the tool developed by V. Salnikov 
et al. in [T]. With this tool using the RSA algorithm generation method, we are able to 
generate three phased RVEs containing up to 100 identical spherical inclusions and up to 
40% in volume fraction. In the cited paper, it has been shown that the method is very 
efficient being able to generate such RVEs in fractions of a second. The described algorithm 
produces a list of inclusions in a ’’vector” form namely for each sphere coordinates of the 
center (in 3D) and the radius. The spheres do not overlap. In the second step, we use a 
voxelization tool ([2]) to obtain a voxelized sample with three phases. We add a coating 
defined by the user for each spherical inclusion. The figure [T] shows a scheme representing 
a section of a sphere with the notations used throughout the article. In the sequel, we note 
I the layer which is the ratio between the coating width and the radius of a sphere r^, the 
layer is thus a number between 0 and 1. On the figures we see one 3D RVE voxelization 
and two of theses sections. 



FIG. 1: Sphere made of a spherical inclusion with a coating parametrized by a layer I 


The setting of the study is the following: we fix the parameters like the volume fraction of 
spheres fgp equal to 30% (which takes into account spherical inclusion and coating) and the 
phases conductivity tensor for each phase equals to kl where A: is a positive number (it takes 
3 different values for the 3 materials and I stands for the identity matrix in dimension 3. 

In this setting, we study how the thermal properties of a composite vary depending on a 
geometrical parameter denominated as the layer keeping in mind that the spherical volume 
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FIG. 2: Example of one RVE meshed with 3 phases generated by RSA method and two sections 
in 128^ voxels representing the top face and one of these center sections 


fraction is fixed here at 30%. For this purpose, we generated samples with different quantity 
of spheres and different resolutions to observe the behaviour of the composite. 

III. COMPUTATIONAL TECHNIQUES 

Consider a representative volume element V with periodic boundary conditions, and 
denote 6{x) the temperature, 0(£c) the heat flux and L(x) the thermal conductivity at any 
point X E V. 

Fourier’s law in the linear case states that: (f){x) = —L{x) V6{x). Moreover, we sup¬ 
pose that we are in stationary phase and without heat source so the heat flux verifies 
div{(f){x)) = 0. Notice that for a composite material, the thermal conductivity tensor 
depends on the point x\ the dependence is governed by microscopic geometry of the sample, 
namely which phase (matrix, inclusion or coating) the point belongs to. 

The volume element V is subjected to an average temperature gradient (V6{x)) = V0 
which induces local temperature gradient V9{x) and heat flux (j){x) inside the RVE. The 
effective constitutive relations of the composite are the relations between the (spatial) 
average flux $ and average temperature gradient V0. 

It is a common practice to introduce a constant reference tensor to solve the problem 
of recovering the local fields and gradients. Introducing this homogeneous linear tensor, we 
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can rewrite the problem as: 


(j){x) = —L^{x)V6{x) +t{x), \/x G V 

div{<p{x)) = 0, 'ix G V (1) 

t(x) = -(L(x) - L^)Ve{x), VxeV 

where r denotes the polarization tensor. 

The solution of (1) can be expressed in real and Fourier spaces, where • denotes the Fourier 
image of • and are the coordinates in the Fourier space, respectively, with the help of 
the periodic Green operator F° associated with we can also write: 

ve{x) = T^{x)Vt{x) VxeV (2) 


or in Fourier space: 

w(^) = r°(^)T(^) ^ 0, w(o) = 0 (3) 


The equations (2) and (3) give the Lippmann-Schwinger integral equations in real and 
Fourier spaces respectively. The Green operator F° is easily expressed and computed in the 
Fourier space by: 



a, 




n 


The Lippmann-Schwinger equation can be solved iteratively using the algorithm based on 
the accelerated scheme proposed by Eyre and Milton in [Hj and written as algorithm 3 in [Ij 
for the elastic case. For the thermal conductivity case, we obtain with adapted notations 
the following algorithm: 
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Algorithm: FFT-based homogenization scheme for thermal conductivity case 
Initialize V6^{x) = V0, fix the eonvergence eriterion aee. 
while (not converged) 

1. Convergence test: 

if (Comp < o-cc) compute <p^{x) = — L{x)V9^(x), (p = FFT{(jF), 

^,, = 1 mm?) turn 

if {ceq < acc) then eonverged, stop 

2. T^(x) = —{L{x) + L^)V6{x) 

3. r” = FFT{t^) 

*■ Cl,JO = -r«(e)f”(?),« 7^ 0, wL,(0) = ve 
5^ = -FFr-‘(V«Lp) 

6. £»™p = ^(||V«” - v«”„,|p)/||ve|| 

7. V9^+\x) = V9^(x) - 2(L(x) - L'^)-^L'^(V9^^p(x) - V9^(x)) 

When this algorithm converges, we compute {<p{x)) and we obtain Lhom using the following 
equality: {(fix)) = -Lhom{'^9{x)) = -Lhom'^S. 

There are several FFT-based numerical schemes like the ’basic scheme’, the ’dual scheme’, 
the augmented Lagrangian scheme, the ’polarization scheme’ and the ’accelerated scheme’. 
We have chosen the last one for its computational efficiency like in the elastic case [Ij. 
We justify this choice by relying on Moulinec and Silva in [H]; they have shown that the 
accelerated scheme is the optimal compromise except for composite with infinite contrast. 
In this paper, the contrast of the materials studied is between 5 and 2.10^. 

IV. COMPARISON FEM AND FFT METHOD 

We briefly describe how we evaluate the homogenized thermal properties of the compos¬ 
ite using the FEM. We use a double-scale method described by Sanchez-Palencia [6j and 
Bensoussan, Lions and Papanicolaou [7]. It deals with setting a multi-scale problem via 
an asymptotic expansion of the equations describing the behaviour of the material. In our 
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case, we consider a double-scale problem in the context of the thermal conductivity. The 
composite material is associated to the macroscopic scale and the RVE to the microscopic 
one. We impose periodic boundary conditions. 

In order to compare both approaches, we study simple cases namely: one sphere or four 
spheres centred as it is shown on hgure 



FIG. 3: Two cuts representing RVE with one coated sphere and four coated spheres in resolution 

64^ voxels with sufficient coating 


We generate RVEs in resolution 64^ and 128^ voxels and with several thickness of the layer. 
For each generated RVE, we construct a mesh based on the voxels then each voxel defines 
a cub8 element and we use a python code to compute the thermal homogenized tensor. 

The RVEs with a thin coating in low resolution namely 64^ voxels or lower, exhibits some 
voids on the coating as one can notice on figure In these cases, the results obtained by 
the two methods show some deviation which can reach 20%. 



L- 


FIG. 4: A cut representing RVE with one sphere with thin coating 


However, for the RVEs with sufficient coating such in figure the results obtained by the 
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different methods are very similar, the discrepancy between both methods does not exceed 
1%. An important remark is that with more complex geometries, the 3D homogenization 
requires a lot of memory resources and takes a lot of CPU time with the FEM whereas the 
EFT method has a lower cost in time and it allows the computation in 256^ voxels and even 
in 512^ voxels (using a computation cluster). Thus we extend the FFT method to more 
complex geometries namely with at least 30 spheres located randomly without intersection. 


V. NUMERICAL SIMULATIONS 

Numerical computations are made for 3 resolutions: 64^, 128^ and 256^ voxels. According 
to these resolutions, for the same geometries, the number of voxels representing each phase 
increases. We also take 30 up to 100 spheres for Hgp. We recall that the layer is a real 
number between 0 and 1. We make the tests with I equal to 0.02 up to 0.08 and to 0.1 up 
to 0.9. In fact, from a layer equal to 0.2, the coating reaches almost the half of the sphere 
volume as we see in the table |T| The coating volume fraction is evaluated by the equality : 
fcoating = (1 ~ (1 “ 0^) ^ where 0.3 denotes the spheres volume fraction. 


For each set of (i?es, Ugp, 1), we generate 10 different samples in order to make the average 
of obtained macroscopic conductivity tensor. 

To compute the thermal conductivity tensor, we recall that the materials are supposed to 
be isotropic namely with a tensor equal to kl, where A: is a positive number and I stands 
for the identity matrix in dimension 3. We take: 


■^matrix phase 


1 1 Lj 


inclusion phase 


= 0 . 2 /, L 


coating phase 


= 400/ 


With theses values, the higher contrast is 2.10^ and the lower is 5. We deliberately choose 
these coefficients for industrial applications with matrix and inclusion phases poorly con¬ 
ducting and coating phase highly conducting. According to [Sj and [S] the constant reference 
thermal tensor is set to: 





min(x) 

€{1,0.2,400} 


X max{x) I 

xe{l,0.2,400} 


-4^/5I 



layer 

fcoating 

0.02 

0.0176 

0.04 

0.0346 

0.06 

0.0508 

0.08 

0.0664 

0.1 

0.0813 

0.2 

0.1464 

0.3 

0.1971 

0.4 

0.2352 

0.5 

0.2625 

0.6 

0.2808 

0.7 

0.2919 

0.8 

0.2976 

0.9 

0.2997 


TABLE I: Coating volume fraction depending on the thickness of the layer 


We remark that the computed thermal homogenized tensor is close to an isotropic mate¬ 
rial tensor, so we take as the homogenized thermal conductivity the mean of the diagonal 
elements of the tensor for the series of the 10 geometries generated as discussed above. 

The curves below on figures ism give these homogenized values. First, we naturally 
observe that the homogenized thermal conductivity increases with the thickness of the layer 
and for all the three resolutions tested. We can also note that the larger the resolution is, 
the less the dispersion is. We conclude that resolution 64 is not accurate enough. 

We observe for very thin layers that the homogenized thermal conductivity decreases when 
the Hgp increases because the coating is thin and the number of voxels representing it decrease 


where we use the equality: 

1 /3 

Number of voxels = I x Res x , with Res G {64,128, 256}. 

The values in the tables are deliberately not rounded up. 


and can be less than one pixel, i.e. it can be missed: this is shown in the tables 
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FIG. 5: Resolution 64^: homogenized thermal conductivity evolution vs layer for risp from 30 up 

to 100 


1 \ 

30 

40 

50 

60 

70 

80 

90 

100 

0.02 

0.17 

0.16 

0.14 

0.14 

0.13 

0.12 

0.12 

0.11 

0.04 

0.34 

0.31 

0.29 

0.27 

0.26 

0.25 

0.24 

0.23 

0.06 

0.51 

0.47 

0.43 

0.41 

0.39 

0.37 

0.36 

0.34 

0.08 

0.68 

0.62 

0.58 

0.54 

0.52 

0.49 

0.47 

0.46 

0.1 

0.86 

0.78 

0.72 

0.68 

0.64 

0.62 

0.59 

0.57 

0.2 

1.71 

1.55 

1.44 

1.36 

1.29 

1.23 

1.19 

1.15 

0.3 

2.57 

2.33 

2.16 

2.04 

1.93 

1.85 

1.78 

1.72 

0.4 

3.42 

3.11 

2.89 

2.72 

2.58 

2.47 

2.37 

2.29 


TABLE II: Voxels coating thickness for the resolution 64^ voxels 
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FIG. 6: Resolution 128^: homogenized thermal conductivity evolution vs layer for Ugp from 30 up 

to 100 


We also note for resolution 256^ voxels that the number of spheres almost does not influence 
the computed values, which agrees with the definition of acceptable RVE size and means we 
no longer have the ’artefact’ due to the voxelization. Moreover, for a layer I equal to 0.08, 
there is a kind of stagnation, it is shown on figures [T^ On these figures, we clearly 
observe two behaviours, one before the stagnation and the other one during the stagnation. 
This effect is less visible in resolution 128 since it is difficult to capture the coating when 
the layer is too low which appears in figures [T^ From the point of view of applications we 
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FIG. 7: Resolution 256^: homogenized thermal conductivity evolution vs layer for risp from 30 up 

to 100 


1 \ 

30 

40 

50 

60 

70 

80 

90 

100 

0.02 

0.34 

0.31 

0.29 

0.27 

0.26 

0.25 

0.24 

0.23 

0.04 

0.68 

0.62 

0.58 

0.54 

0.52 

0.49 

0.47 

0.46 

0.06 

1.03 

0.93 

0.87 

0.81 

0.77 

0.74 

0.71 

0.69 

0.08 

1.37 

1.24 

1.15 

1.09 

1.03 

0.99 

0.95 

0.92 

0.1 

1.71 

1.55 

1.44 

1.36 

1.29 

1.23 

1.19 

1.15 

0.2 

3.42 

3.11 

2.89 

2.72 

2.58 

2.47 

2.37 

2.29 

0.3 

5.13 

4.66 

4.33 

4.07 

3.87 

3.70 

3.56 

3.44 


TABLE III: Voxels coating thickness for the resolution 128^ voxels 
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1 \ 

30 

40 

50 

60 

70 

80 

90 

100 

0.02 

0.68 

0.62 

0.58 

0.54 

0.52 

0.49 

0.47 

0.46 

0.04 

1.37 

1.24 

1.15 

1.09 

1.03 

0.99 

0.95 

0.92 

0.06 

2.05 

1.87 

1.73 

1.63 

1.55 

1.48 

1.42 

1.37 

0.08 

2.74 

2.49 

2.31 

2.17 

2.06 

1.97 

1.90 

1.83 

0.1 

3.42 

3.11 

2.89 

2.72 

2.58 

2.47 

2.37 

2.29 


TABLE IV: Voxels coating thickness for the resolution 256^ voxels 


can say that it is not necessary to have a big coating to increase the homogenized thermal 
conductivity. There is a limit due to the sphere volume fraction and due to the type of 
spherical inclusions. 


2.8 

2.6 

2.4 
2.2 
2.0 
1.8 
1.6 

1.4 
1.2 
1.0 
0.8 
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FIG. 8: Homogenized thermal conductivity evolution vs coating volume fraction for 30 spheres 


VI. CONCLUSIONS AND OUTLOOK 

To conclude, we can summarize the main ideas of this paper. We have shown that 
the FFT-based iterative accelerated scheme is a good tool for computing thermal effective 
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FIG. 9: Homogenized thermal conductivity evolution vs coating volume fraction for 70 spheres 
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FIG. 10: Homogenized thermal conductivity evolution vs coating volume fraction for 100 spheres 


properties of a composite generated randomly in 3D. We have restricted our study to RVE 
composed with only spheres but it can be made with cylinders or both spheres and cylinders. 
In [1], RVEs generation algorithms in MD (Molecular dynamics) or RSA methods have 
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FIG. 11: Example of sections with lack of coating when the layer is too low 


been presented nsing both spheres and cylinders and in jl] these RVEs have been nsed for 
calculations in elastic case. 

We have observed that the resolution 256^ with 50 — 70 spheres and an sufficient thickness 
of the layer gives good results. We have also noted that the coating must be thick enough 
in order to be ’’captured”. In our work in progress, we explore the case with a thin coating, 
where the methods are more subtle and some preparation of the samples is needed. We 
will address this issue as well as more advanced morphological parameters in a separate 
publication. 
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